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A technique for estimating aerodynamic parameters in real time from flight data without 
air flow angle measurements is described and demonstrated. The method is applied to 
simulated F-16 data, and to flight data from a subscale jet transport aircraft. Modeling 
results obtained with the new approach using flight data without air flow angle 
measurements were compared to modeling results computed conventionally using flight data 
that included air flow angle measurements. Comparisons demonstrated that the new 
technique can provide accurate aerodynamic modeling results without air flow angle 
measurements, which are often difficult and expensive to obtain. Implications for efficient 
flight testing and flight safety are discussed. 
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body-axis translational accelerometer measurements, g 
wing span, ft 

mean aerodynamic chord, ft 

body-axis nondimensional aerodynamic force coefficients 
body-axis nondimensional aerodynamic moment coefficients 
expectation operator 
Fourier transform operator 
mass moments of inertia, slug-ft 2 

imaginary number = V-T 
cost function 
aircraft mass, slug 

body-axis pitching moment from engine thrust, ft-lbf 
body-axis roll, pitch, and yaw rates, rad/s or deg/s 
ambient static pressure, lbf/ft 2 
dynamic pressure, lbf/ft 2 

specific gas constant for air = 1716.59 ft-lbf/slug-deg R 

wing reference area, ft 2 

maneuver length, s 

ambient temperature, deg R 

body-axis engine thrust, lbf 

body-axis air-relative velocity components, ft/s 

airspeed, ft/s 

angle of attack, rad or deg 
sideslip angle, rad or deg 

stabilator, elevator, aileron, and rudder deflections, rad or deg 
Euler roll, pitch, and yaw angles, rad or deg 
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= parameter vector 
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p 

= density, slug/ft 3 

a 

= standard error 

X 

= covariance matrix 

CO 

= frequency, rad/s 

superscripts 

T 

= transpose 

f 

= complex conjugate transpose 

* 

= estimate 


= time derivative 

~ 

= Fourier transform 

-1 

= matrix inverse 

subscripts 

o 

= reference value 


I. Introduction 

One practical challenge associated with estimating aerodynamic model parameters from flight data is obtaining 
accurate measurements of the air flow angles, namely the angle of attack and sideslip angle. Instrumenting an 
aircraft for the air flow angles is not easy, because the introduction of any sensor modifies the local air flow, and the 
local air flow is also modified by the presence of the aircraft itself. Air flow angle sensors are often moved away 
from the aircraft using wingtip booms or a nose boom, to reduce local flow effects. Because flight data analysis for 
aircraft stability and control is based on air flow angles at the aircraft center of gravity, air flow angle measurements 
from sensors mounted on booms must be corrected to the center of gravity using the measured rotational rates of the 
aircraft and the position of the sensor relative to the aircraft center of gravity 1 . Hypersonic vehicles cannot use air 
flow angle sensors that protrude into the air flow, because the aerodynamic heating that attends hypersonic flight 
would destroy the sensors. On many commercial transport aircraft, air flow angle sensors are mounted on the 
fuselage. This location requires extensive and costly calibration for good accuracy. 

Flight vehicle aerodynamics depend very strongly on angle of attack and sideslip angle, so it might at first seem 
impossible to conduct flight data analysis and modeling without direct measurements of these quantities. In the past, 
short-term data compatibility analysis 1 and various filtering techniques 2 have been used to reconstruct air flow 
angles from other measurements. However, these methods induce drift in the reconstructed quantities over time, as 
a result of integrating kinematic equations with inputs that are aircraft motion measurements containing bias errors. 

Many aircraft can be inexpensively equipped with inertial measurement units located fairly close to the aircraft 
center of gravity. Instrumentation for air flow angles is much more difficult and expensive, both to install and to 
calibrate properly. It would be advantageous if aerodynamic model parameters could be determined accurately from 
flight data without having to instrument the aircraft to measure air flow angles. 

In this paper, a method is proposed which is based on the general principle of reconstructing air flow angles from 
inertial data, while taking advantage of the characteristics inherent in a frequency-domain parameter estimation 
technique, to obtain accurate aerodynamic model parameter estimates without direct measurements of the air flow 
angles. Estimated aerodynamic model parameters include static stability derivatives associated with angle of attack 
and sideslip angle, along with all other traditional stability and control derivatives. The approach is an extension of 
a real-time parameter estimation technique 1,3 ' 4 that has been successfully demonstrated in flight 5 ' 7 , and also 
independently evaluated 8,9 as an excellent method for real-time dynamic modeling. The new approach can be 
implemented in real time, or as a post-flight method applied to all the data at once. 

Although the proposed method has obvious advantages for reducing costs associated with stability and control 
flight testing, there are important implications for aircraft safety. Commercial and general aviation aircraft typically 
do not have calibrated, research-quality air flow angle sensors mounted on booms, but many do have 3-axis inertial 
measurement units or could be equipped with that instrumentation inexpensively. The proposed method extends 
real-time aircraft dynamic modeling capability to these aircraft, which makes the idea of an onboard, real-time 
stability and control advisory to the pilot a realistic prospect. Quantitative information about the stability and 
control of the aircraft, provided to the pilot in real time, would improve aircraft safety and could also be used for 
adaptive or reconfigurable control, or monitoring aircraft degradation due to aging or ordinary wear. With the 
proposed method, commercial aircraft could have this capability with just an update to non-flight-critical software, 
or with added hardware that could passively access flight data that already exists on the onboard data bus. 
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Section II describes the method. The model formulation retains full nonlinear aircraft dynamics, with linearized 
aerodynamic models. Parameters in the linear aerodynamic models can be time varying, because the modeling 
algorithm provides frequent updates to the model parameter estimates in real time. Descriptions of the F-16 
nonlinear simulation and the subscale jet transport aircraft are given in Section III, including flight instrumentation 
and characteristics of the flight data. The test aircraft has air flow angle sensors mounted on booms attached to each 
wingtip. This makes it possible to do the aerodynamic parameter estimation in a conventional way using air flow 
angle measurements, and compare results with aerodynamic parameter estimates from the proposed method, which 
does not use air flow angle measurements. Flight results are presented in Section IV. Finally, the concluding 
remarks in Section V summarize the technical contribution and applications. 

II. Methods 


A. Air Flow Angle Reconstruction 

The translational equations of motion for a rigid aircraft in body axes are 1 : 


qSC x . T 

u = rv — qw -l — g sin 0 + — 

(la) 

m m 

qSCy „ . , 

v = pw -ru-i — + g cos 0 sin <p 

(lb) 
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qS 

w = qu - pv -t —+ g cosdcostp + — 

(lc) 


m m 

Substituting translational acceleration measurements for the applied forces results in the translational kinematic 
equations in body axes 1 : 


it = rv-qw— g sinO + ga x 

(2a) 

v = pw - ru + geos 0 sin <f>+ ga v 

(2b) 

w = qu — pv + g cos 6 cos f + ga z 

Airspeed, sideslip angle, and angle of attack can be computed from u, v, and w using 1 

(2c) 

1/ / 2,2,2 
V =sj u + v +w 

(3a) 

P = sin '(v/V) 

(3b) 

a = tan '(w/m) 

(3c) 


Many flight test maneuvers for aerodynamic parameter estimation are small perturbation maneuvers about a 
reference condition of steady, straight-and-level flight at low nominal angle of attack, in low- wind conditions. With 
these assumptions, the following approximations are valid: 

u ~V ~ constant (5 « v/V a « w/V (4) 

For subsonic flight, airspeed can be computed from measurements of dynamic pressure and ambient conditions, 

V =^2q/p 
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P = Pal{RT a ) 


so that 

V=pqRT a /P a 

Combining Eqs. (2b), (2c), and (4), 


P ~ pa- r + ^(cosd sin</> + a y J 


a&q-pP+ ^(cosOcostfr + a.) 


For a trim condition of steady, symmetric, straight-and-level flight. 


u = v = P = w = a = 0 
p = q= r = 0 
v = P = a y =0 = O 
0 = a 


( 5 ) 


(6a) 

(6b) 


( 7 ) 


Combining these trim conditions with Eqs. (2), initial conditions for /? and a at the start of a maneuver are 

P{ 0) = 0 (8a) 

a(0) = sin 1 ( a Xg j (8b) 

Eqs. (5)-(8) can be used to reconstruct sideslip angle /? and angle of attack a from other measurements. The 
practical problem with this approach is that the measured values for the inertial data, namely the angular rates 
p,q, and r , and the translational accelerometer measurements a x , a v , and a z , all typically have bias errors. This 

causes a time-dependent drift in the P and a reconstructions from Eqs. (6), even when the bias errors are small, 
because of the additive effect of the time integration. Any bias error in the steady-state accelerometer measurement 
a Xo or inaccuracy in the level trim condition will bias the initial condition for the angle of attack reconstruction, and 

any deviation from symmetric flight in the initial trim condition means that using zero as the initial condition for 
sideslip angle reconstruction will be incorrect. Scale factor errors for angular rate and translational accelerometer 
measurements are negligible in practice, and any errors in the Euler angles <f> and 0 will be mitigated by the 
trigonometric functions and the division by airspeed, which is a relatively large value for an airplane. Many inertial 
measurement units provide Euler angle time histories (or equivalent quaternion data) that have been corrected for 
drift, using magnetometer measurements, for example. Failing that, the Euler angles can be reconstructed from 
angular rate data by integrating the rotational kinematic equations 1 

(j> = p + tan 9 [q sin <j> + r cos (f) (9a) 

0 = q cos (f>- r sin (j) (9b) 


with initial conditions 
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< P { 0)=0 


(10a) 


&(0) = a(0) = sin l (a x j 


(10b) 


From the foregoing discussion, it is clear that reconstructed J3 and a time series can have significant bias and 
drift errors. Flowever, if the modeling is done in the frequency domain, where the bias and drift are removed for 
other reasons (discussed next), then this method of reconstructing air flow angle data can be used very effectively. 

When transforming time-domain data into the frequency domain, the constant bias and drift (also called the 
trend, or the best-fit linear function of time) are always removed prior to applying the Fourier transform. This 
processing is done to avoid spillage from relatively large-amplitude low-frequency components that can pollute the 
frequency-domain data at low frequencies of interest 1 ' 4 . It follows that bias and drift errors in the air flow angle 
reconstructions do not have to be corrected when the analysis is done in the frequency domain, because the bias and 
drift of every time series are thrown out anyway. In essence, the reconstruction is done in the time domain, and the 
errors incurred in doing that are discarded prior to transformation into the frequency domain. At that point, the 
analysis proceeds as usual, using frequency-domain transform methods. Note that this approach takes advantage of 
the fact that angular rate measurements and translational accelerometer measurements typically have small but 
non- zero bias errors, and negligible slope or scale factor errors 1 . 

Although the focus here is on airplane problems, a similar approach could be used for rotorcraft, by using 
Eqs.(2) in place of Eqs. (6) for the data reconstruction. This could be particularly useful in hover, where body-axis 
velocity components u, v, and w are relatively small, and body-axis velocity components generally cannot be 
measured accurately because of local flow effects from the rotor downwash. 

The next two subsections explain how the aerodynamic modeling can be done in real time in the frequency 
domain. Similar material can be found in Refs. [l],[3]-[7]. 


B. Aerodynamic Modeling 

Nondimensional aerodynamic force and moment coefficients for an aircraft can be computed from flight 
measurements as follows 1 : 


C = 


)-m 7 


C x =(ma x -T x )/qS 
C z = (ma z -T z )/q s 

M + {I x - ] z ) P r + f xz {P 2 ~ r2 ) 

C Y = ma y /qS 

Q = [ r xP - r xz {pq + r) + ( ■ I z - Jy ) qr ] /qSb 

c„ = [ I z i '- I xz(p-q r )+{iy- I x)pq]/qSb 


qSc 


(11a) 

(lib) 

(11c) 

(12a) 

(12b) 

(12c) 


These expressions retain the frill nonlinear dynamics in the equations of motion for a symmetric rigid aircraft. 
For local modeling over a short time period, the force and moment coefficients computed from Eqs. (11) and (12) 
can be modeled using linear expansions in the aircraft states and controls: 


Cx 

Cz 


— C x A a C x 

= Cy Aa+Cy 

Z. a z. 


Aqc 
a 2V 

Aqc 
« 2V 


- C x AS + C x 


- Cy AS + Cy 


(13a) 

(13b) 
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(13c) 


Aqc 


C m — C w Aa + C... 1- C m AS + C m 

m m a m q 2 y mg m ° 


A pb Arb 

Cy — Cy Afl + Cy + Cy + Cy AS + Cy 

1 l p *p 2V Ir 2V 15 Iq 

(14a) 

C, = C t Aj3+ C, Aph +C, Arb + C, AS + C, 

1 ip ip 2 Y l r 2 Y 15 l ° 

(14b) 

C n =C n A/3 + C n Apb + C n Arb +C n A S + C„ 

n np r n p ^ y n r ^ y n s n 0 

(14c) 


The A notation indicates perturbation from a reference condition. In each expansion, a single term is shown to 
represent all relevant and similar control terms, to simplify the expressions. For example, in Eq. (14b), the term 
C, AS represents all the control terms for C, , e.g., Q AS = C, A8 a + C, A8 r . In Eq. (14b), Q represents 

O S Sq 5y o 

the nondimensional rolling moment at a reference condition, and similarly for the other expansions. 

The linear aerodynamic models in Eqs. (13) and (14) contain parameters called stability and control derivatives, 
such as Q /; and Cj g , that characterize the stability and control of the aircraft. For short periods of time, the 

stability and control derivatives are considered as constant model parameters to be estimated from flight data. 

The next subsection describes how the unknown stability and control derivatives in the linear aerodynamic 
models of Eqs. (13) and (14) can be estimated from flight data using equation-error parameter estimation in the 
frequency domain. 

C. Stability and Control Derivative Estimation in the Frequency Domain 

The first step required for modeling in the frequency domain is to transform the measured flight data from the 
time domain into the frequency domain. The finite Fourier transform is the analytical tool used for this task. For an 
arbitrary scalar signal x(t ) on the time interval [0, 7’] , the finite Fourier transform is defined by 

lT[A'(f)] = T(o) = J o x(t) e~^ co, dt (15) 


which can be approximated by 

N - 1 

x(®) « At ^ x(i) e~^ mAt (16) 

i '=0 


where x(i) = x(iAt) , T = (lV-l)zlf , and At is a constant sampling interval. The summation in Eq. (16) is 
defined as the discrete Fourier transform, 


v-i 

X(®)= Yj x (0 e ~ j03iA> ( 17a ) 

i '=0 

so that the finite Fourier transform approximation in Eq. (16) can be written as 

x(®) « X (co) At (17b) 

Some fairly straightforward corrections 10 can be made to remove the inaccuracy resulting from the fact that 
Eq. (17b) is a simple Euler approximation to the finite Fourier transform in Eq. (15). Flowever, if the sampling rate 
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is much higher than the frequencies of interest, as is typically the case for dynamic modeling from flight data, then 
the corrections are small and can be safely ignored. 

The Fourier transform is applied to the nondimensional force and moment coefficients computed from Eqs. (11) 
and (12) using measured time-domain data. This results in the nondimensional force and moment coefficients in the 
frequency domain. Often, measurements of the angular accelerations p, q, and r are not available. In the frequency 
domain, these derivatives can be calculated by multiplying the Fourier transforms of p,q, and r by jco . For 
example, the Fourier transform of the rolling moment coefficient from Eq. (12b) can be computed as: 


C/H 


T[C,\ = 


jcoT 


I X P - I X7. r 

+ 3= 

-txzPq + 1 

jz-'y) 

1 qr 

qSb 

qSb 



(18) 


and similarly for C m and C n . This approach implements the derivative of the body-axis angular momentum in the 
frequency domain, including the time variation in the inertia quantities. Note that the Fourier transform of the 
nonlinear terms is handled by computing the nonlinear terms in the time domain, then applying the Fourier 
transform to the resulting time history. Treatment of the dynamic pressure q in Eq. (18) is consistent with an 
assumption that the dynamic pressure varies slowly, which is a good practical assumption. 

To obtain the perturbation states and controls in Eqs. (13) and (14), time histories of the measured states and 
controls are high-pass filtered to remove the bias (steady part) and drift (linear function of time) for each signal. 
Eqs. (5)-(8) are used to reconstruct sideslip angle and angle of attack from other measurements. The time series for 
reconstructed sideslip angle and angle of attack are high-pass filtered in the same manner as the other data, resulting 
in perturbation time series for sideslip angle and angle of attack. All perturbation signals are then transformed into 
the frequency domain using the discrete Fourier transform. The break frequency for the high-pass filter is set below 
the lowest frequency selected for the modeling. Similarly, the quantities transformed in Eq. (18) (shown within the 
square brackets) are high-pass filtered prior to Fourier transformation. This approach effectively drops out the bias 
terms in the models of Eqs. (13) and (14), and allows the use of a simple multiplication by jco for differentiation in 
the frequency domain. The high-pass filtering also preserves the accuracy of spectral data near zero frequency by 
avoiding leakage from the relatively large signal components associated with bias and drift. For reconstructed 
/:> and a , the low-frequency components removed by the high-pass filtering are bias and drift errors resulting from 
the time-domain reconstruction. 

When a time series has non-zero endpoints, the Fourier transform of the time derivative of that time series 
requires endpoint correction terms. Applying integration by parts to Eq. (15), 


= J x(t)e ■' c0> dt 

= x(t) e~j coT - x(jS) + jcjj x(t)e ^'"dt (19a) 

= x(T ) e ^" T - x(0) + jco x(a>) 

Applying high-pass filtering to a time series imposes a zero initial condition. This leaves only one term for the 
endpoint corrections, so that 


^'[x(t)] = x(t) e ■ 7<ur + jcox(co) (19b) 

Eq. (19b) was used to compute the finite Fourier transform for the time derivative of a time series that has been 
high-pass filtered. 

For each aerodynamic model in Eqs. (13) and (14), the parameter estimation problem can be formulated in the 
frequency domain as a standard least squares regression problem with complex data 1 , 


Z = XO + e 


(20) 
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where, for example, using the rolling moment equation (14b), 


Q(i) 

C,( 2 ) 

C,(M) 




A Pn{ 1 ) 

A r n { 1 ) 

A *a{ 1 ) 

A S r (l) ' 

x = 

Af(2) 

A Pn( 2 ) 

^,( 2 ) 

^ 4 ( 2 ) 

A *r(2) 


A~P(M ) 

A Pn(M) 

A U M ) 

A S a {M) 

A5 r {M ) 


0 = 


C, 


c, 

l p 

c, 


c, 


c, 


(21) 


(22) 


(23) 


The notation Ap n represents T [ pb / 2V ] , and similarly for Ar n , and e represents the complex vector of equation 
errors in the frequency domain. The symbols Afi(k), k = 1,2, M denote the Fourier transform of the sideslip 
angle perturbation state for each frequency co k , and similarly for other quantities. Each transformed variable 
depends on frequency. The frequencies o:> k can be chosen arbitrarily, and are therefore chosen to cover the 
frequency band where the aircraft dynamics lie, as will be discussed later. The least squares cost function is 


J{e) = ±(i-xe)\z-xe) 


(24) 


This cost function contains M squared-error terms in summation, corresponding to M frequencies of interest. 
Similar cost expressions can be written for individual lines from Eqs. (13) and (14). The parameter vector estimate 
that minimizes the least squares cost function is computed from 1 



X f X 




(25) 


The estimated parameter covariance matrix is 1 





2 


= a 



where the equation-error variance a 2 can be estimated from the residuals. 
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(27) 


V — 

( M ~ n p 

and n p is the number of unknown parameters, i.e., the number of elements in parameter vector 9 . Estimated 

parameter standard errors are computed as the square root of the diagonal elements of the X ( 0 ) matrix from 

Eq. (26), using d 2 from Eq. (27). Explanations of why the estimated parameter standard errors are computed in 
this way, and why this calculation in the frequency domain produces parameter error measures that are consistent 
with the scatter in parameter estimates from repeated maneuvers, can be found in Ref. [1], Realistic simulation 
testing and flight test data analysis have shown that the accuracy of model parameters estimated with this method in 
real time is comparable to using a time-domain output-error method employing iterative nonlinear optimization in 
post-flight batch mode 1 ' 4 ' 11 . 

The model formulation given here is widely applicable, because the assumption of constant linear aerodynamic 
models over short time periods is very accurate for nondimensional stability and control derivatives, where the 
effects of changing dynamic pressure and mass properties are removed. Nonlinear model terms, such as 

C m ^Aa 2 or C m Aa Aqc /2V , could also be used in the models of Eqs. (13) and (14) without any change to the 

a 2 aq 

procedure described here. The nonlinear terms would simply be computed in the time domain, using perturbation 
signals generated from measured signals by high-pass filtering to remove the bias and drift. The resulting nonlinear 
model term would then be transformed into the frequency domain, like any other time series. Therefore, equation- 
error parameter estimation in the frequency domain can also be used with nonlinear model terms, although the 
models would still be linear in the unknown parameters. 

To implement this least squares parameter estimation in the frequency domain, the parameter estimation 
calculations in Eqs. (25)-(27) are applied to frequency-domain data at selected times, normally at regular time 
intervals, such as 0.5 s, which corresponds to an update rate of 2 Hz. Other values for the time interval could also be 
chosen, but linearized aerodynamic characteristics can normally be captured very well by updating the parameter 
estimates at a rate of 1 or 2 Hz, except in cases of strong nonlinearity, damage, failure, or rapid maneuvering. For 
those situations, the update rate can be increased, at the cost of more computations. The frequency-domain data 
must therefore be available at any time, so the Fourier transforms are computed using a recursive Fourier transform, 
described next. 

D. Recursive Fourier Transform 

For a given frequency co , the discrete Fourier transform in Eq. (17a) at time i At , denoted by X t (to) , is related 
to the discrete Fourier transform at time (i - l)At by 


) 


(z-xdf [z-xe] 


Xi(a>) = X i _ l (a>) + : x(i)e ]mAt 


where 


e ~ja>iAt _ e -j<oAt e ~jco(i-\)At 


(28) 


(29) 


The quantity e ja>Al is constant for a given frequency co and constant sampling interval At . It follows that the 
discrete Fourier transform can be computed for a given frequency at each time step using one addition in Eq. (28) 

and two multiplications - one in Eq. (29) using the stored constant e~ ja>Al for frequency co , and one in Eq. (28). 
There is no need to store the time-domain data in memory when computing the discrete Fourier transform in this 
way, because the data for each sample time is processed immediately. Time-domain data from the past can be used 
in all subsequent analysis by simply continuing the recursive calculation of the Fourier transform. In this sense, the 
recursive Fourier transform acts as memory for the information in the data. More data from more maneuvers 
improves the quality of the data in the frequency domain without increasing memory requirements to store it. 
Furthermore, the Fourier transform is available at any time i At . The approximation to the finite Fourier transform 
is completed using Eq. (17b). 
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The recursive computation of the Fourier transform does not use a Fast Fourier Transform [FFT ) algorithm 12 , 
and therefore would be comparatively slow, if the entire frequency band up to the Nyquist frequency 1 /2 At were of 
interest. Flowever, rigid-body dynamics of full-scale aircraft lie in a rather narrow frequency band of approximately 
0.01-2.0 FIz. The frequency range shifts higher for subscale aircraft 1 , but the width of the frequency band is similar. 
Since the frequency band is limited, it is efficient to compute the discrete Fourier transform using Eqs. (28) and (29) 
(which represents a recursive formulation of Eq. (17a)) for only the selected frequencies a\, k =1,2 . With 
this approach, it is possible to select closely-spaced fixed frequencies for the Fourier transform and subsequent 
modeling, and still do the calculations efficiently. 

Using a limited frequency band for the Fourier transformation confines the data analysis to the frequency band 
where the dynamics lie, and automatically filters wide band measurement noise or structural responses outside the 
frequency band of interest. This is important for real-time applications, because both transformation into the 
frequency domain and automatic noise filtering are accomplished with the simple calculations in Eqs. (28) and (29). 

In past work on fighter aircraft short-period modeling, frequency spacing of 0.04 FIz on an interval of 
approximately [0.1-2] FIz was found to be adequate 3 ' 7 . Finer frequency spacing requires slightly more computation, 
but was found to have little effect on the results. When the frequency spacing is very coarse, there is a danger of 
omitting important frequency components, and this can lead to inaccurate parameter estimates. In general, a good 
rule of thumb is to use frequencies evenly spaced at 0.04 Hz over the bandwidth for the dynamic system. For good 
results, the bandwidth should be limited to the frequency range where the signal components in the frequency 
domain are at least twice the amplitude of the wide band noise level. However, the algorithm is robust to these 
design choices, so the selections to be made are not difficult. 

The recursive Fourier transform update need not be done for every sampled time point. Systematically skipping 
time points effectively lowers the sampling rate of the data prior to Fourier transformation. This saves computation, 
and, assuming proper analog anti-alias filtering 1 has been implemented in the instrumentation system, does not have 
a significant adverse impact on the parameter estimation results until the Fourier transform update rate gets below 
approximately 5 times the highest frequency of interest for the dynamic system. The parameter estimation and 
covariance calculations in Eqs. (25)-(27) can be done at any time, but are usually done at 1 or 2 Hz, to save 
computations. Linearized aerodynamic characteristics rarely change faster than this, except in cases of strong 
nonlinearity, damage, failure, or rapid maneuvering. For these cases, the update rate can be increased, at the cost of 
more computations. 

Estimated parameter standard errors computed from the covariance matrix in the frequency domain using 
Eqs. (26) and (27) do not require correction for colored residuals 1 . These standard errors are therefore a good 
representation of the uncertainty in the estimated parameters. 

III. Aircraft 

A. F-16 Nonlinear Simulation 

The F-16 is a single-seat, multi-role fighter with a blended 
wing / body and a cropped delta wing planform with leading edge 
sweep of 40 deg. Thrust is provided by one General Electric FI 10- 
GE-100 or Pratt & Whitney F100-PW-220 afterburning turbofan 
engine mounted in the rear fuselage. Figure 1 shows a photograph 
of the F-16. Mass properties and geometry of the aircraft are given 
in Table 1. 

The aircraft was modeled with controls for throttle s th > 
stabilator 8 S , aileron 8 a , and rudder 8 r . Speed brake and flaps 
were assumed fixed at zero deflection. Throttle deflection was 
limited to the range 0 <8 th < 1 , stabilator deflection was limited to 

-25° < 8 S < 25° , aileron deflection was limited to 

-21.5° < 8 a <21.5° , and rudder deflection was limited to -30° <8 r < 30° . These limits represent the physicals 
stops. 



Figure 1. F-16 Aircraft 
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Non-dimensional aerodynamic force and moment coefficient data were derived from a low-speed static wind 
tunnel test and a dynamic forced oscillation wind tunnel test, both conducted on a 16% scale model of the F-16. The 
aerodynamic database applies to the F-16 flown out of ground effect, with landing gear retracted, and no external 
stores. Static aerodynamic data are in tabular form as a function of angle of attack and sideslip angle over the ranges 

-10° < a < 45° and -30° < /? < 30° , respectively. Dynamic data is provided in tabular form at zero sideslip angle 

over the angle of attack range -10° < a < 45° . Dependence of the non-dimensional coefficients on a is included 
in the q dependencies, due to the manner in which the data was collected in the wind tunnel. 

The F-16 nonlinear simulation is programmed completely in MATLAB ’ . Full nonlinear equations of motion, 
including turbine engine gyroscopic effects, are used. Complete details on the F-16 nonlinear simulation can be 
found in Appendix D of Ref. [1]. 


B. T-2 Subscale Jet Transport Aircraft 

The T-2 aircraft is a 5.5 percent subscale model of a 
generic commercial twin-engine jet transport aircraft. A 
photograph of the aircraft in flight is shown in Figure 2. 

The aircraft has twin jet engines mounted under the 
wings and retractable tricycle landing gear. Mass 
properties and geometry of the aircraft are given in 
Table 1. Further information on the T-2 subscale jet 
aircraft and associated flight test operations can be found 
in Refs. [13] and [14], 

1 . Control Surfaces 

Control surfaces on the aircraft are ailerons, inboard and outboard elevators, upper and lower rudders, inboard 
and outboard trailing-edge flaps, and inboard and outboard wing strakes, for a total of 16 independent control 
surfaces. For the data analyzed in this work, only the elevators, ailerons, and rudders were deflected. The 
individual elevator surfaces were moved together as a single elevator surface, and similarly for the rudders. Left and 
right ailerons were deflected asymmetrically, in the conventional way. Definitions of control surface deflections are 
given below. Trailing edge down is positive deflection for the wing and elevator surfaces, and trailing edge left is 
positive for the rudder. 



Figure 2. T-2 Subscale Jet Transport Aircraft 
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The aircraft can be flown by a safety pilot using direct visual contact and conventional radio control. A research 
pilot executes flight test maneuvers from inside a mobile control room, using a synthetic vision display drawn from 
telemetry data and a local terrain database. Inputs from the research pilot and a ground-based flight control system 
are used to compute control surface commands which are transmitted by telemetry to the aircraft. 

2. Instrumentation and Data Acquisition 

The T-2 aircraft was equipped with a micro-INS, which provided 3-axis translational accelerometer 
measurements, angular rate measurements, estimated attitude angles, and GPS velocity and position. Air data 
probes attached to booms mounted on each wingtip (visible in Figure 2) measured angle of attack, sideslip angle, 
static pressure, and dynamic pressure. Measurements from static pressure sensors and ambient temperature sensors 
were used to compute air density and altitude. Engine speeds in rpm were measured and used as inputs to an engine 
model to compute thrust. The engine model was identified from ground test data, with adjustments for ram drag 
identified from flight data. Potentiometers on the rotation axes of the control surfaces measured control surface 
deflections. Mass properties were computed based on measured fuel flow, pre-flight weight and balance, and inertia 
measurements of the aircraft on the ground. The pilot stick and rudder pedal commands and throttle position were 
also measured and recorded. Data from onboard sensors were telemetered to the ground in real time. Sampling rate 
for the flight data was 50 FIz. 
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IV. Results 


A. Simulation Results 

The method for aerodynamic parameter estimation using reconstructed air flow angle data was first applied to 
data from a nonlinear simulation of the F-16 aircraft 1 . A doublet sequence on the stabilator, rudder, and ailerons 
was applied to the simulation, and Eqs. (5)-(8) were used to reconstruct air flow angle data using simulated 
measurements from other sensors. Five percent Gaussian random noise was added to all simulation outputs, and 
typical bias errors (approximately 0.1 deg/s for the angular rate sensors and 0.01 g for the translational 
accelerometers) were applied. Figure 3 shows the comparison of simulated measured sideslip angle and angle of 
attack from the maneuver and the reconstructions using other measurements. The time-domain reconstructions 
clearly show bias and drift errors, but the frequency domain comparison is so close that the traces are almost 
indistinguishable. This occurs because high-pass filtering in the time domain effectively removed the bias and drift 
in both the measured and reconstructed air flow angle data prior to Fourier transformation. Some experimentation 
and investigation into the characteristics of high-pass filters indicated that good results could be obtained using a 
4 th -order Butterworth high-pass filter, with cut-off frequency set at f mjn / 2 , where f min is the minimum frequency 
used in the frequency-domain modeling. Other high-pass filtering solutions could likely be devised to work as well. 

Figure 4 shows that frequency-domain estimates for the model parameters associated with a and /tare accurate 
even when reconstructed data for sideslip angle and angle of attack are used. The symbols marked as “simulation” 
indicate parameter estimates computed using finite differences applied to the aerodynamic database in the F-16 
nonlinear simulation, with perturbation sizes chosen to match the variable ranges exhibited during the maneuver. 
The symbols marked “measured” are the parameter estimates using direct measurements of sideslip angle and angle 
of attack from the simulation, and the symbols marked “reconstructed” are parameter estimates based on sideslip 
angle and angle of attack data reconstructed from other measurements using Eqs. (5)-(8). In both cases, parameter 
estimation results were computed using the real-time frequency-domain method described in Section II. The small 
vertical bars through the “measured” and “reconstructed” symbols represent 95 percent confidence intervals (±2cr) 

for the parameter estimates, computed from Eqs. (26) and (27). The results show that the real-time frequency- 
domain method produced accurate parameter estimates from reconstructed air flow angle data, i.e., without air flow 
angle measurements. Estimates of stability and control derivatives not associated with a and f3 were similarly 
unaffected by whether or not reconstructed air flow angle data were used. 

The F-16 nonlinear simulation generated simulated air flow angle measurements at the aircraft center of gravity. 
In practice, measurements from air flow angle sensors located away from the aircraft center of gravity would have to 
be corrected to the aircraft center of gravity prior to stability and control modeling. Similarly, accelerometer 
measurements must be corrected to the aircraft center of gravity. When using reconstructed air flow angle data in 
the frequency domain, only the accelerometer position correction is needed, because air flow angle measurements 
are not used. This results in a computational savings that partially offsets the computations required for air flow 
angle data reconstruction. 

A Monte Carlo analysis was done by selecting angular rate and translational acceleration biases from specified 
intervals using a uniform probability density, and applying white Gaussian noise sequences to the simulation 
outputs. The angular rate biases were chosen on the interval [-0.01, 0.01] rad/s and translational acceleration biases 
were chosen on the interval [-1.0, 1.0] ft/s 2 . Parameter estimation results were computed using air flow angle data 
from the F-16 simulation, and also using reconstructed air flow angle data, as in the single run used to generate 
Figures 3 and 4. The real-time frequency-domain method described in Section II was used to do the modeling for 
each Monte Carlo run, and the entire analysis was repeated 200 times with the same doublet input sequence applied 
to the control surfaces, but with different instrumentation biases and random noise sequences for each run. 
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Figure 3. F-16 Simulated Data Reconstruction 
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Figure 4. Static Stability Derivative Estimates from F-16 Simulation Data 


Table 2 shows a comparison of the parameter estimation results from the Monte Carlo runs. The standard errors 
were computed from the scatter in the parameter estimates, 


a g = 



(31a) 


(31b) 


where n is the total number of Monte Carlo runs. The computed a g quantifies the scatter in repeated parameter 
estimates. 

Results in Table 2 show that the frequency-domain modeling approach using reconstructed ah flow angle data is 
accurate in terms of the mean value of the parameter estimates, but generally exhibits greater scatter in the 
individual estimates from each run, compared to using the same parameter estimation method and measured air flow 
angle data. However, the magnitudes of the scatter in the parameter estimates using reconstructed air flow angle 
data are still relatively small. Figure 5 shows histograms of the Monte Carlo parameter estimation results for one 

particular model parameter, C z . The figure shows that the mean value of C z using reconstructed air flow angle 
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data agreed with the mean value of C z using measured air flow angle data, but the scatter increased using 

reconstructed air flow angle data. In both cases, the scatter in repeated estimates was relatively small, with 2a 
values less than 2 percent of the mean value. Estimates for other model parameters showed similar behavior. 

Larger scatter in the parameter estimates using reconstructed air flow angle data makes sense from a data 
information perspective, because the missing air flow angle sensor data means that less information is being 
supplied to the estimator. In any case, error bounds for the parameter estimates can be lowered by running 
maneuvers repeatedly and averaging the individual parameter estimates. Because the proposed method is accurate 
in the mean, improved parameter accuracy can be obtained without measurements of the air flow angles by simply 
using repeated maneuvers to compensate for the loss of information from the missing air flow angle sensor 
measurements. 

Parameter estimation results from the Monte Carlo analysis demonstrate that accurate modeling results can be 
obtained using the proposed method with reconstructed air flow angle data and realistic bias errors in the angular 
rate and translational acceleration measurements. 






Figure 5. F-16 C 7 Derivative Estimates for 200 Monte Carlo Simulation Runs 


All parameter estimation results (i.e., in Figures 3-5 and Table 2) were computed using the real-time formulation 
of equation-error parameter estimation in the frequency domain, as described in Section II. The parameter 
estimation results shown are the final values from the real-time estimator, available at the end of the maneuver. The 
behavior of real-time parameter estimates and error bounds at other times during the maneuver was investigated in 
Ref. [7] using flight data and various inputs. 


B. Flight Results 

After the method was successfully developed and validated using F-16 nonlinear simulation data, the same 
approach was applied to flight data from the subscale T-2 aircraft. This section contains the flight results and 
associated discussion. 

Figure 6 shows flight data from the T-2 aircraft during a maneuver where 3-axis optimized multi-sine 
inputs 1 ' 7 ' 15 ' 16 were applied simultaneously to the elevator, aileron, and rudder controls. The nominal flight condition 
was 5 deg angle of attack and 130 ft/s airspeed, at 1250 ft altitude. Figure 7 shows a comparison between static 
stability derivative estimates computed using sideslip angle and angle of attack measurements from vanes on the 
wingtip booms, along with estimates of the same parameters computed using sideslip angle and angle of attack data 
reconstructed from Eqs. (5)-(8). In both cases, real-time parameter estimation in the frequency domain, as explained 
in Section II, was used to produce the modeling results. The results shown in Figure 7 are the final values from the 
real-time parameter estimator, at the end of the maneuver. These values match batch estimates computed using the 
same frequency-domain least-squares method, which is equivalent to stating that the recursive calculation of the 
discrete Fourier transform in Eqs. (28)-(29) is an exact equivalent (within the numerical precision of the 
computations) to the batch computation given in Eq. (17a). The symbols in Figure 7 indicate the estimated 
parameter values, and the vertical lines represent 95 percent confidence intervals (±2<x). Estimates for the other 
stability and control derivatives showed consistency similar to that shown in Figure 7. 
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Figure 6. T-2 Flight Data for 3-axis Optimized Multi-Sine Inputs 
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Figure 7. Static Stability Derivative Estimates from T-2 Flight Data 
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As in the F-16 simulation case, the parameter estimates computed using reconstructed air flow angle data 
showed good consistency but slightly larger uncertainties, compared to the parameter estimates computed using 
measured air flow angle data. Figure 8 shows a comparison of high -pass filtered time series for measured and 
reconstructed air flow angles. Reconstruction of the air flow angle data was not perfect, but captured the main 
frequency content, which is most important. The mismatches in the reconstructions shown in Figure 8 include the 
effects of imperfections in the air flow angle measurements from the wingtip vanes, such as the structural response 
of the wings and wingtip booms during dynamic maneuvering, turbulence, inertia and friction in the vane 
mechanisms, and local flow effects, typically called sidewash and upwash. The modeling results shown here 
demonstrate some robustness to an imperfect match between measured and reconstructed air flow angle data. This 
type of realistic mismatch was not modeled in the F-16 simulation case. 

An important measure of the quality of an identified model is prediction capability for data not used in the 
modeling process. Figure 9 shows a prediction case where the dynamic model identified using reconstructed air 
flow angle data was used to predict flight data from a 3-axis doublet maneuver at a similar flight condition. 
Although the plots indicate a slight error in predicting the Dutch roll natural frequency, the overall prediction 
accuracy is very good. Because the parameter estimates using measured air flow angle data agreed well with the 
parameter estimates using reconstructed air flow angle data, the slight prediction error seen in Figure 9 also occurs 
when using a model with parameters estimated using measured air flow angle data. 
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Figure 8. Sideslip Angle and Angle of Attack Reconstructions from T-2 Flight Data 
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Figure 9. T-2 Flight Data Prediction for 3-Axis Doublet Inputs 
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V. Concluding Remarks 

Simulated data and flight data were used to demonstrate that real-time aerodynamic parameter estimation can be 
done in the frequency domain without measurements of air flow angles, i.e., sideslip angle and angle of attack. Air 
flow angle data was reconstructed in the time domain, then detrended using high-pass filtering, and transformed into 
the frequency domain for modeling. Necessary assumptions were that flight test maneuvers must involve small 
amplitude perturbations, so that a simple form of the air flow angle data reconstruction equations can be used, and 
that the angular rate and translational acceleration data have bias errors but not scale factor errors. Both of these 
conditions are typical in practice. Once the air flow angle data reconstruction is done, the data are high-pass filtered, 
then transformed into the frequency domain, and the data analysis and modeling proceeds as described in previous 
works, including the capability to accurately estimate aerodynamic model parameters and error bounds in real time. 

The work presented here had its genesis in necessity, from hypersonic flight tests and subscale aircraft flight 
tests with limited instrumentation. In these situations, direct measurement of air flow angles was either not 
practical, or expensive in terms of both time and money. The method described here can be particularly useful for 
subscale model testing with limited resources, where the installation of an inertial measurement unit close to the 
center of gravity is relatively inexpensive and easy to do. Adding a pitot tube for static and dynamic pressure, and 
instrumentation to measure control surface deflections, results in an aircraft instrumented for useful data analysis 
and dynamic modeling in terms of dimensional stability and control derivatives, using the frequency-domain method 
described in this paper. To estimate nondimensional stability and control derivatives, mass properties of the aircraft 
must also be determined, typically from ground tests or from detailed computer-aided design (CAD) modeling of the 
aircraft mass distribution. 

Results from simulation investigations showed that the technique described here could be used effectively for 
accurate real-time estimation of aircraft stability and control parameters. Validation was done by comparing results 
computed using the proposed method without air flow angle data to results from a conventional aerodynamic 
parameter estimation approach using the air flow angle measurements, and to values obtained from finite differences 
applied to the simulation aerodynamic database. Results from a subscale jet transport aircraft showed that the 
proposed real-time frequency-domain method works very well in practice, and gives aerodynamic parameter 
estimates comparable in accuracy to those obtained using air flow angle measurements from vane sensors installed 
on wingtip booms. The dynamic models identified without air flow angle measurements showed excellent 
prediction capability for a maneuver with dissimilar input forms at a similar flight condition. Modeling results can 
be obtained in real time, using a high-pass filter and the recursive Fourier transform, as described in the paper. 

The parameter estimation approach described here is simple and accurate, and produces high-quality modeling 
results for flight testing with limited instrumentation. The method extends the capability to provide high-quality 
real-time stability and control information to aircraft that have poor measurements of air flow angles, or do not have 
air flow angle instrumentation at all. This capability has important implications for improved aircraft safety, 
because the method could be implemented in any aircraft that has a flight computer and data from an inertial 
measurement unit, along with ambient temperature and pressure, dynamic pressure, and measured control surface 
deflections. Many modern commercial aircraft have this data available at sample rates of 25 Hz, meaning that 
stability and control monitoring could be provided to the pilot in real time with only a software change to non-flight- 
critical software in the flight control computer, or by adding hardware that accesses flight data available from the 
onboard data bus. 
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Table 1. F-16 and T-2 Geometry and Mass Properties 



F-16 

T-2 

length c , ft 

11.32 

0.908 

wing span b , ft 

30 

7.083 

wing area S , ft 2 

300 

7.046 

x ref , in 

0.35 c 

42.628 

y ref ’ ^ 

0 

0.000 

z re f i in 

0 

0.000 

x cg , in 

0.25 c 

42.728 

y cg > in 

0 

0.000 

z C g > in 

0 

0.519 

m , slugs 

637.16 

1.502 

I x , slugs-ft 2 

9,496 

1.077 

ly , SlllgS-ft 2 

55,814 

4.163 

1 z , slugs-ft 2 

63,100 

5.016 

I xz , slugs-ft 2 

982 

0.416 
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Table 2. Monte Carlo Parameter Estimation Results for the F-16 Nonlinear Simulation 


for 200 Maneuvers using Sequential Doublet Inputs 
V 0 = 450 ft/sec, h B = 10,000 ft, a Q = 0 O = 5 deg 


Parameter 

e 

Simulation 
F inite- Difference 
Value 

Measured a, fi 
Real-Time 
Frequency-Domain 
Estimate O + a 

Reconstructed a, fi 
Real-Time 
Frequency-Domain 
Estimate 0 + a 

C x 

0.3101 

0.3161 ±0.0013 

0.3 149 ±0.0022 

Cx « 

0.0000 

0.6444 ±0.0540 

0.6515 ±0.1099 

Cx S s 

0.0812 

0.0607 ±0.0020 

0.0616 ±0.0034 

C z 

-3.6154 

-3.6104 ±0.0127 

-3.5952 ±0.0251 

c z 

-31.400 

-31.334 ±0.519 

-31.473 ± 1.138 

c z 

-0.4354 

-0.4414 ±0.0194 

-0.4511 ± 0.0378 

r 

m 

" l a 

-0.3443 

-0.3476 ±0.0036 

-0.3463 ±0.0041 

r 

m q 

-8.4000 

-8.1129 ±0.1298 

-8.1204 ±0.1478 

r 

m S s 

-0.5926 

-0.5788 ±0.0058 

-0.5797 ±0.0065 

Cy 

Y P 

-1.1459 

-1.1443 ±0.0051 

-1.1358 ± 0.0539 

Cy 

1 P 

0.1100 

0.1175 ±0.0295 

0.1367 ±0.2079 

Cy 

± r 

0.9580 

0.9673 ±0.0625 

0.9846 ±0.1121 

Cy 

Y Sa 

0.0602 

0.0629 ±0.0085 

0.0688 ±0.0686 

Cy 

Y S r 

0.1642 

0.1639 ±0.0024 

0.1635 ±0.0044 

c, 

V 

-0.1375 

-0.1352 ±0.0007 

-0.1342 ±0.0063 

Cl 

l p 

-0.4200 

-0.4087 ±0.0038 

-0.4065 ± 0.0245 

Cl 

l r 

0.1130 

0.1180 ±0.0073 

0.1200 ±0.0136 

c, 

c a 

-0.1490 

-0.1450 ±0.0013 

-0.1443 ±0.0082 

C 1 

l s r 

0.0267 

0.0256 ±0.0004 

0.0255 ± 0.0006 

c n 

n p 

0.2610 

0.2594 ±0.0014 

0.2577 ±0.0123 

C n 

n p 

-0.0162 

-0.0206 ±0.0081 

-0.0241 ±0.0473 

c n 

n r 

-0.4421 

-0.4167 ±0.0176 

-0.4196 ±0.0246 

C n 

fix 

°a 

-0.0281 

-0.0289 ± 0.0027 

-0.0300 ±0.0155 

C ”s 

°r 

-0.0921 

-0.0896 ±0.0009 

-0.0896 ±0.0012 
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